!****************************************************************************
! Below are routines to approximately translate between R and nu, where
!
! nu := [delta_c(z)/(G(z)*sigma(R))]**2,
!
! delta_c(z) is the threshold density for spherical collapse (Henry 2001)
! G(z) is what I am calling the growth factor (e.g. Peebles 1980 eq. 10.12)
! sigma**2(R) is the variance of the linear power spectrum top-hat smoothed
! on a scale R.
!
! nu0(m) := (delta_c,0/sigma(m))**2
!
! delta_c,0 := 3(12PI)**(2/3)/20 ~= 1.68, from Henry (2001), is the
!      constant part of the spherical collapse density threshold. So
!      delta_c(z)=delta_c,0*(1-0.0123*LOG((1/OMEGA_M0-1)a**3). The extra
!      term depending on z and cosmology is added whenever
!      nu=(delta_c(z)/(G(z)*sigma(R))**2 is calculated. The function
!      Delta_cFactor(z) (below) calculates 1-0.0123*LOG((1/OMEGA_M0-1)a**3.
!      I am calling the growing factor G(z).
!
! sigma**2(m) is the variance of the linear power spectrum top-hat
!      smoothed over a scale R(m)=[3*m/(4*PI*avg_density)]**(1/3).
!      It is normalized.
!
!****************************************************************************
module mNu

        use mPrecision
        use mObject
        use mLPS
        use mFit

        implicit none

        interface new
                module procedure new_nu
        end interface new

        interface getUniverse
                module procedure getUniverse_nu
        end interface getUniverse

        interface getLPS
                module procedure getLPS_nu
        end interface getLPS

        interface nu_m
                module procedure nu_m_nu
        end interface nu_m

        interface m_nu
                module procedure m_nu_nu
        end interface m_nu

        interface evaluate
                module procedure nuEvaluate
                !module procedure evalChebyshevFunction
        end interface evaluate

        interface evaluateInverse
                module procedure nuEvaluateInverse
        end interface evaluateInverse

        interface fitLL
                module procedure fitLL_nu
        end interface fitLL

        interface fitUL
                module procedure fitUL_nu
        end interface fitUL

        interface dnu
                module procedure dnu_nu
        end interface dnu

        real(D), parameter :: LLM = 0.9900e3, &
                              ULM = 1.1e16

        private
        public :: new
        public :: getUniverse
        public :: getLPS
        public :: m_nu, nu_m
        public :: dnu
        public :: Delta_cFactor
        public :: refit

contains

        subroutine new_nu(nu, lps, nu0ADependence, a)

                type(tNu), intent(out) :: nu
                type(tLPS), intent(in), optional :: lps
                logical, intent(in), optional :: nu0ADependence
                real(D), intent(in), optional :: a

                if( present(lps)) then
                        nu%lps = lps
                else
                        call new(nu%lps)
                end if

                if( present(nu0ADependence) ) then
                        nu%nu0ADependence = nu0ADependence
                else if( transferFunction(nu%lps) == 'EH' ) then
                        nu%nu0ADependence = .true.
                else
                        nu%nu0ADependence = .false.
                end if

                if( present(a) ) then
                        nu%a=a
                else
                        nu%a=1.0_D
                end if

                call new(nu%nuFit, LLM, ULM, inverse = .true., &
                         xLog = .true., yLog = .true., &
                         param = nu%a, lps = nu%lps)
                call fit(nu%nuFit, nu0)

                nu%dnuFit_a=-1.0_D

        end subroutine new_nu

        !********************************************************************
        !   nu0(M) := (delta_c,0/sigma(M))**2
        !
        !   delta_c,0 := 3(12PI)**(2/3)/20 ~= 1.68, from Henry (2001), is the
        !      constant part of the spherical collapse density threshold. So
        !      delta_c(z)=delta_c,0*(1-0.0123*LOG((1/OMEGA_M0-1)a**3). The
        !      extra term depending on z and cosmology is added whenever
        !      nu=(delta_c(z)/(G(z)*sigma(R))**2 is calculated. The function
        !      Delta_cFactor(z) (below) calculates
        !      1-0.0123*LOG((1/OMEGA_M0-1)a**3. I am calling the growing
        !      factor G(z).
        !
        !   sigma**2(M) := the variance of the linear power spectrum top-hat
        !      smoothed over a scale R(M)=[3*M/(4*PI*avgDensity)]**(1/3)
        !
        !*******************************************************************
        real(D) function nu0(params, M)

                use mConstants
                use mUniverse
                use mLPS
                use mPSFittingFunctions

                type(tFitParams), intent(in) :: params
                real(D), intent(in) :: M
                real(D) :: a

                a = param(params)
#ifdef NON_LINEAR_MASS_function
                nu0 = 9.0_D * (12.0_D * PI) ** (4.0_D / 3.0_D) &
                      / (400.0_D * smithSmoothed( getLPS(params), &
                                    m_R(universe( getLPS(params) ), M), a ))
#else ! NON_LINEAR_MASS_function
                nu0 = 9.0_D * (12.0_D * PI) ** (4.0_D / 3.0_D) &
                      / (400.0_D * smoothedVariance( getLPS(params), &
!                                    m_R(universe( getLPS(params) ), M), a, 'TH'))
                                    m_R( getLPSUniverse(params), M ), a, 'TH'))
#endif ! NON_LINEAR_MASS_function
        end function nu0

        type(tUniverse) pure function getUniverse_nu(nu)
                type(tNu), intent(in) :: nu
                getUniverse_nu = getUniverse( getLPS(nu) )
        end function getUniverse_nu

        type(tLPS) pure function getLPS_nu(nu)
                type(tNu), intent(in) :: nu
                getLPS_nu = nu%lps
        end function getLPS_nu

        real(D) pure function fitLL_nu(nu)
                type(tNu), intent(in) :: nu
                fitLL_nu = rFitLL(nu%nuFit)
        end function fitLL_nu

        real(D) pure function fitUL_nu(nu)
                type(tNu), intent(in) :: nu
                fitUL_nu = rFitUL(nu%nuFit)
        end function fitUL_nu

        real(D) pure function nu0LL(nu)
                type(tNu), intent(in) :: nu
                nu0LL = evaluate( nu, fitLL(nu) )
        end function nu0LL

        real(D) pure function nu0UL(nu)
                type(tNu), intent(in) :: nu
                nu0UL = evaluate( nu, fitUL(nu) )
        end function nu0UL

        real(D) pure function nuEvaluate(nu,x)
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: x
                nuEvaluate = evalRationalFunction(nu%nuFit, x)
        end function nuEvaluate

        real(D) pure function nuEvaluateInverse(nu,x)
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: x
                nuEvaluateInverse = fitEvaluateInverse(nu%nuFit, x)
        end function nuEvaluateInverse

        !**********************************************************************
        !* Approximates normalized nu := [delta_c(z) / (G(z) * sigma(R))] ** 2.
        !**********************************************************************
        real(D) function R_nu(nu,R,a)
                use mUniverse
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: R, a

                call checkA(nu, a)
                                
                in_domain: if(R < m_R( getUniverse(nu), fitLL(nu) ) &
                              .or. R > m_R( getUniverse(nu), fitUL(nu) )) then
                        write(0, *) 'ERROR: R_nu(R,a) called with R =', R, &
                                    'but this is out of range:'
                        write(0, *) '             ', &
                                m_R( getUniverse(nu), fitLL(nu) ), &
                                '<= R <=', &
                                m_R( getUniverse(nu),fitUL(nu) )
                        stop
                end if in_domain
                R_nu = nuADependence(nu, a) ! DEREK---square this
                R_nu = R_nu * R_nu * evaluate(nu, R_m( getUniverse(nu), R ))
        end function R_nu

        !********************************************
        ! Approximates R given a normalized
        ! nu := [delta_c(z) / (G(z) * sigma(R))] ** 2
        !********************************************
        real(D) function nu_R(nu, nuValue, a)
                use mUniverse
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: nuValue, a
                real(D) :: nu0

                call checkA(nu, a)
                nu_R = nuADependence(nu, a)
                nu_R = nu_R * nu_R
                nu0 = nuValue / nu_R
                in_domain: if( nu0 < nu0LL(nu) .or. nu0 > nu0UL(nu) ) then
                        write(0, *) 'ERROR: nu_R(nu,a) called with nu=', &
                                     nuValue, ' but this is out of range:'
                        write(0, *) '     ', nu0LL(nu) / nu_R, '<= nu <=', &
                                     nu0UL(nu) / nu_R
                        stop
                end if in_domain
                nu_R = m_R( getUniverse(nu), evaluateInverse(nu, nu0) )
        end function nu_R

        !**************************************************************
        ! Approximates normalized
        ! nu := [delta_c(z) / (G(z) * sigma(R(m)))] ** 2.
        !
        ! R(m) := (3 * m / (4 * PI * rho^bar)) ** (1 / 3) is calculated
        ! by m_R from module mMisc
        !**************************************************************
        real(D) function m_nu_nu(nu, m, a)
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: m, a
                call checkA(nu, a)

                in_domain: if( m < fitLL(nu) .or. m > fitUL(nu) ) then
                        write(0, *) 'ERROR: m_nu(mass,a) called with mass=', &
                                     m ,'but this is'
                        write(0, *) '       out of range:', fitLL(nu), &
                                    '<= mass <=', fitUL(nu)
                        stop
                end if in_domain

                m_nu_nu = nuADependence(nu, a)
                m_nu_nu = m_nu_nu * m_nu_nu * evaluate(nu, m)
        end function m_nu_nu

        !****************************************************
        ! Approximates normalized m given
        ! nu := [delta_c(z) / (G(z) * sigma(R(m)))] ** 2.
        ! m(R) := 4 / 3 * PI * rho^bar * R ** 3 is calculated
        ! by R_m from module mMisc
        !****************************************************
        real(D) function nu_m_nu(nu,nuValue,a)
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: nuValue, a
                real(D) :: nu0

                call checkA(nu, a)
                nu_m_nu = nuADependence(nu, a)
                nu_m_nu = nu_m_nu * nu_m_nu
                nu0 = nuValue / nu_m_nu

                in_domain: if( nu0 < nu0LL(nu) .or. nu0 > nu0UL(nu) ) then
                        write(0, *) 'ERROR: nu_m(nu,a) called with nu=', &
                                     nuValue
                        write(0, *) '       but this is out of range'
                        write(0, *) '     ', nu0LL(nu) * nu_m_nu, &
                                    '<= nu <=', nu0UL(nu) * nu_m_nu
                        stop
                end if in_domain

                nu_m_nu = evaluateInverse(nu, nu0)
        end function nu_m_nu

        !************************************************************
        ! Access to raw fitting function, for informational purposes.
        !************************************************************
        real(D) function R_nu0(nu,R)

                use mUniverse

                type(tNu), intent(in) :: nu
                real(D), intent(in) :: R
                in_domain: if(R < m_R( getUniverse(nu), fitLL(nu) ) &
                    .or. R > m_R( getUniverse(nu), fitUL(nu) )) then
                        write(0, *) 'ERROR: R_nu0(R,a) called with R=', R, &
                                    '       but this is out of range:'
                        write(0, *) '      ', m_R( getUniverse(nu), fitLL(nu) ), &
                                    '<= R <=', m_R( getUniverse(nu), fitUL(nu) )
                        stop
                end if in_domain
                R_nu0 = evaluate(nu, R_m( getUniverse(nu), R ))
        end function R_nu0

        real(D) function nuADependence(nu, a)
                use mUniverse
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: a
                nuADependence = Delta_cFactor( getUniverse(nu), a ) &
                    / growthFactor( getUniverse(nu), a )
        end function nuADependence

        subroutine checkA(nu, a)
                type(tNu), intent(in) :: nu
                real(D), intent(in) :: a
                if(nu%nu0ADependence) then

                        if(a /= nu%a) then
                            write(0, *) 'WARNING: nu fit function being used'
                            write(0, *) '         for the wrong scale factor'
                        end if

                end if
        end subroutine checkA

        subroutine refit(nu,a)
                type(tNu), intent(inout) :: nu
                real(D), intent(in) :: a
                if(nu%nu0ADependence) then
                        if(a /= nu%a) then
                                write(0, *) 'INFO: Refitting nu'
!                                call destroyNu(nu)
                                call new(nu, nu%lps, a = a)
                                nu%nu0ADependence = .true.
                        end if
                end if
        end subroutine refit

        real(D) function dnu_nu(nu, M, a)
                type(tNu), intent(inout) :: nu
                real(D), intent(in) :: M, a
                call checkDNuFit(nu, a)
                dnu_nu = m_nu(nu, M, a) / M * log( evaluate(nu%dNuFit, M) )
        end function dnu_nu

        subroutine checkDNuFit(nu, a)
                type(tNu), intent(inout) :: nu
                real(D), intent(in) :: a
                if(nu%nu0ADependence) then
                        if(a /= nu%a .or. a /= nu%dNuFit_a) then
                            write(0, *) 'WARNING: nu fit function being used'
                            write(0, *) '         for the wrong scale factor'
                        end if
                end if
                if(nu%dNuFit_a == -1.0_D) call initDNuFit(nu, a)
        end subroutine checkDNuFit

        subroutine initDNuFit(nu, a)
                type(tNu), intent(inout) :: nu
                real(D), intent(in) :: a
                type(tCFit) :: tempFit
                write(0, *) 'Fitting dnu'
                call new(tempFit, LLM, ULM, xLog = .true., yLog = .true., &
                         param = a, lps = nu%lps)
                call fit(tempFit, nu0)
                nu%dNuFit = chebyshevDerivative(tempFit)
                nu%dNuFit_a = a
        end subroutine initDNuFit

        !**********************************************************************
        ! Calculates the redshift/cosmology dependent part of delta_c, the
        ! spherical collapse threshold, using Henry's (2001) fitting formula.
        ! I have a separate function for this because the constants are already
        ! contained in my rational function approximation for nu.
        !**********************************************************************
        real(D) pure function Delta_cFactor(u, a)
                use mUniverse
                type(tUniverse), intent(in) :: u
                real(D), intent(in) :: a ! scale factor a=1/1+z
                Delta_cFactor = 1.0_D - 0.0123_D &
                    * log(1.0_D + ( 1.0_D / getOmega_m0(u) - 1.0_D ) * a ** 3)
        end function Delta_cFactor

end module mNu
